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The Rayleigh-Taylor instability of two immiscible fluids in the limit of small Atwood 
numbers is studied by means of a phase-field description. In this method the sharp fluid 
interface is replaced by a thin, yet finite, transition layer where the interfacial forces vary 
smoothly. This is achieved by introducing an order parameter (the phase field) whose 
variation is continuous across the interfacial layers and is uniform in the bulk region. The 
phase field model obeys a Cahn-Hilliard equation and is two-way coupled to the standard 
Navier-Stokes equations. Starting from this system of equations we have first performed 
a linear analysis from which we have analytically rederived the known gravity-capillary 
dispersion relation in the limit of vanishing mixing energy density and capillary width. 
We have performed numerical simulations and identified a region of parameters in which 
the known properties of the linear phase (both stable and unstable) are reproduced in 
a very accurate way. This has been done both in the case of negligible viscosity and in 
the case of nonzero viscosity. In the latter situation only upper and lower bounds for 
the perturbation growth-rate are known. Finally, we have also investigated the weakly- 
nonlinear stage of the perturbation evolution and identified a regime characterized by a 
constant terminal velocity of bubbles/spikes. The measured value of the terminal velocity 
is in perfect agreement with available theoretical prediction. The phase-field approach 
thus appears to be a valuable tecnhique for the dynamical description of the stages where 
hydrodynamic turbulence and wave-turbulence enter into play. 



1. Introduction 

The Rayleigh-Taylor (RT) instability is a fluid-mixing mechanism occurring when a 
heavy, denser, fluid is pushed into a lighter one. For a fluid in a gr avitational field , 



such a mechanism was first discovered by Lord Rayleigh in the 1880s (Rayleigh 
and later applied to all accelerated fluids by Sir Geoffrey Taylor in 1950 (|Tavlor 
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The relevance of this mixing mechanism embraces many different phenomena occur- 
ring in completely different contexts. We just mention, among the many, astrophysi- 
cal su pernova explosions and geophysical formations like salt domes and volcanic is 



cal su pernova explosion s and ge ophysical tor mations like salt domes and volcanic is- 
lands (|Di Prima fc Swinnev lll98ll: lDimonte fc Schneider Il2000h. continental magmatism 



caused by lithospheric gravitational instability (|Lee. Rudnick fc Brimhall Jr. l200ll : lDucea fc Saleebv 
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Il998h . inertial c onfinement fusion (Cook fc Zhou 20021 ) and cloud formation in atmo- 



spheric sciences l|Schultz et al. 112009 1 



Back to classical fluids applications, RT instability is the first step eventually leading to 
a fully developed turbulent regime. A deeper understanding of the mechanism of flows 
driven by RT instability thus would shed light on the many processes that underpin fully 
developed turbulence. 

The difficulty inherent in sustaining an unstable density stratification has challenged 
experimentalists for over half a century. Several innov ative approaches have been recently 



developed (see e.g., iRamaprabhu fc Andrews 1120041) 



With the advent of supercomputers, high-resolution numerical simulations of RT at high 
Reynolds numbers have become a reality. However, simulations using many different 
benchmark codes and experiments disagree already on apparently innocent observables 
like, for instance, the value of t he growth constant, a. assoc iated to the spread of the 
turbulent mixing zone (see, e.g., |Pi Prima fc Swinnev 1119811) . The differences can be as 
high as 100%. 

Despite the long history of RT tur bulence, a c onsist ent phenomenological theory has 
been presented only very recently by IChertkov I ([2003) for the miscible case. Th e theo- 



retical predictions by Chertkov have been verified by Celani. Mazzino fc Vozella ( 20061 ) 



exploiting numerical simulati ons in two spatial d imensions. For the three-dimensional 



miscible case we refer, e.g, to lYoung et al. ( 200l[ ) 



In many of the aforementioned situations where the RT instability has an impor- 
tant role, the two fluids are immiscible owing to a non negligible surface tension. At 
level of linear analysis the role played by a non zero surface tension was addressed by 



Chandrasekhar ( 196 if). The successive dynamics falling in a turbulent regime has been 

^ (|2005l) . Using a phenomenological 



recently analyzed bv lChertkov. Kolokolov fc Lebedev 



approach, the authors suggest the existence of a Kolmogorov cascade between the integral 
scale and a time-dependent scale related to the typical drop size. Below the latter scale, 
associated to an emulsion-like region, a wave energy cascade takes in. This is mediated 
by weakly interacting capillary waves propagating on top of the drop surface. Eventually, 
the energy is dissipated by viscous forces. 

RT instability and RT turbulence of immiscible fluids thus appear richer than the cor- 
responding miscible situations. The existence of two different cascades poses a serious 
challenge to numerical investigations of the immiscible RT problems. The emulsion-like 
phase indeed occurs at very small scales and the energy transfer takes place on the inter- 
faces. These are geometrical objects close to singularities and thus difficult to describe 
appropriately in a numerical scheme. Accuracy and efficiency are thus fundamental re- 
quirements to reproduce the correct statistical features characterizing immiscible RT 
turbulence. 

Our aim here is to perform a first step along this direction by focusing on direct nu- 
merical simulations of immi scible RT instability. The numerical strategy we exploit here is 

known as phase-field mode l (jBrav 2002t lCahn fc Hilliard Hl958l ; lBadalassi. Ceniceros fc Baneriee 
2003t iDing. Spelt fc Shull2007l) . The main idea of the method is to treat the interface 
between two immiscible fluids as a thin mixing layer across which physical properties 
vary steeply but continuously. The evolution of the mixing laye r is ruled by an order pa - 



rameter (the phase field) that obeys a Cahn-Hilliard equation (jCahn fc Hilliard 1 119581 ). 



The method permits to avoid a direct tracking of the interface and easily produces the 
correct interfacial tension from the mixing layer free energy. 

We present here an accurate numerical study that validates the phase-field approach 
by testing known results of immiscible RT instability both at level of linear and weakly 
nonlinear analysis. From our results, it turns out that this strategy is a valuable option for 
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Figure 1. Fluids configuration corresponding to a heavier fluid of density p2 placed above a 

lighter one of density p\ < p2- 

a quantitative treatment of the turbulent regime characterized by the interplay between 
hydrodynamic and interface degrees of freedom. 

The paper is organized as follows. In Seed] we introduce the Rayleigh-Taylor problem 
and discuss the related phase-field approach. A detailed analysis of the energy balance 
between purely hydrodynamic degrees of freedom and interface degrees of freedom is 
presented. Finally, the dispersion relation for gravity-capillary waves is obtained by an- 
alytical calculations starting from the phase-field equation coupled to the Navier-Stokes 
equations. 

In Sec. [3] the results from the direct numerical simulations are presented and compared 
with known results for the linear analysis. We focus both on the case of zero viscosity and 
on that of negligible viscosity. Both stable and unstable configurations are considered. 
Finally, the weakly nonlinear regime is considered and the resulting terminal velocity of 
bubbles/spikes compared with existing theoretical predictions. 
Sec. 2] is devoted to some conclusions and perspectives. 



2. System configuration and phase-field model 

Our system consists of two immiscible, incompressible fluids (labeled by 1 and 2) having 
different densities, p\ and p 2 > pi, with the denser fluid placed, e.g., above the less dense 
one (see Fig. [TJ. In the absence of gravity, this flow configuration is stable. In presence 
of the gravitational force, surface tension may be able to keep the system in equilibrium, 
provided the density contrast is not too large. 

Let us start by describe the equilibrium configuration and then pass to the evolution (RT 
instability) that occurs when a perturbation is imposed to the interface separating the 
two fluids. 

2.1. Equilibrium state 

Let us consider an equilibrium state where fluid 1 is placed below fluid 2 and they 
are separated by a sharp interface. The fact that the interface is sharp (i.e. a discon- 
tinuity in the fluid properties) poses a serious challenge to numerical simulations. In- 
deed, for sharp interfaces, the evolution equations are obtained by following fluid 1 and 
2 separately with the appropriate boundary cond ition at interface (see, for instance, 
Smolianski. Haario fc Luukka 112005c ISethian I ll999h . Other approaches follow the inter- 



face alone. In this latter case, the movement of the interface is naturally amenable to a 
Lagrangian description, while the bulk flow is conventionally solved in an Eulerian frame- 
work. These approaches employ a mesh that has grid points on the interfaces and deforms 
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according to the flow. A major shortcoming of these approaches is in that they cannot 
handle properly topological changes such as breakup, coalescence and reconnections (see 
Yue et al. 112004 . and references therein). In this respect, the phase-field method is, by 
far and large, more effective, at the expense of a larger number of grid points required. 

The idea of the phase-field method is to replace the sharp interface with a diffuse one in 
such a way that the numerical computation of interface movement and deformation ca n 
be carried out on fixed grids (jAnderson. McFadden fc Wheeler 1 1 19981: IJacamin I [l999h . 
More quantitatively, this amounts to assigning to the sy stem a Ginzburg-Landau free 



energy, espressed in term of the order parameter <fi as (jCahn fc Hilliard Nl958t iBrav 



2002; lYue et al. 1 12004') 



(2.1) 



where f2 is the region of space occupied by the system, A is a mixing energy density and 
e is the capillary width, representative of the interface thickness. The order parameter cj> 
is a field which serves to identify fluid 1 and 2. We assume (j> — 1 in the region occupied 
by fluid 1 and <f> = — 1 in those where fluid 2 is present. 

The equilibrium state is the minimizer of the free energy T . The mechanism which keeps 
the system in this configuration is the competition between two effects due to the two 
addends in ()2. 1|) . The first term favours a perfect mixing (i.e. A\d(j>\ 2 /2 = in T, this 
term being the interface energy contribution) whereas the second one one drives the 
system towards demixing (the associated term in J 7 , the bulk contribution, has indeed 
a minimum for <f> = ±1). The nontrivial final equilibrium state is just the results of this 
competition. More quantitatively, the final state is obtained by minimizing the free-energy 
functional with respect to variations of the function 0, i.e., solving: 



[i = 5T/5<i> = & -d 2 (j) + 







(2.2) 



where ix is the so-called chemical potential (see, for instance, Cahn fc Hilliard ll 1958 : Brav 
20021: lYue et al. 1 12004) . If one considers an one-dimensional interface, varying along the 
gravitational direction y, one easil y finds the solution of Eq. (|2.2[) as (jCahn fc Hilliard 
1958HBravl[200l : lYue et al. 1120041) : 



4>{y) = ±tanh ^-J=- 



(2.3) 



This solution exists and is stable in all d imensions although the decay rate of pert ur- 



20051) . From 



bations depends upon the dimensionality (jKorvola. Kupiainen fc Taskinen 
(|2.2[) one immediately realizes that the sharp-interface limit is obtained for e — > 0: in this 
case tanh (y/ (\/2e)) — > sign(y). Moreover, the surface tension a is equal to the integral o f 



the free-energy density along the interfa ce (see, for examp le, Landau fc Lifshitz I 2000) 



For a plane interface, this integral yields ( Cahn fc Hilliard 19581 : IBrav I 20021 : Yue et al 
20041) : 



2a/2 A 



(2.4) 



It is now easy to verify how the sharp interface limit is obtained: it suffices to take 
the limits A and e to zero keeping a fixed to the value prescribed by surface tension 
dLiu fc Shenll2008h . 
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2.2. Perturbation evolution 

Let us now suppose to impose a small perturbation on the (finite thickness) interface 
separating the two fluids. Such perturbation will displace the phase field from the previous 
equilibrium configuration, which minimized the free-energy T , to a new configuration for 
which in general, p ^ 0. The system will react so as to try to reach again an equilibrium 
configuration. In formulae: 



dt 



dej) = 7 a> = j Ad 2 



(2.5) 



7 being the so-called mobility (see, for instance, Brav 12002 ; Yue et al. 2004 ). Notice 
the presence of the Laplacian operator in front of p. Notice that the mass of each fluid 
is conserved, as imposed by the physics of the problem under consideration. 
The dyna mics of the velocity field is governed by the usual Boussinesq Navier-Stokcs 
equations (jKundu fc Cohen 1 l200lh plus an additional stre ss contribut i on arising at the 
interface where th e effect of surface tension enters into play ijBrav 2002; Yue et al. 1 120041 : 



Berti et al. 1120051 ). The equations of motion are: 



{d t v 



v ■ dv C] 



9aP 
Po 



vd 2 v a 



-d, 



8F 



Po 6cj) 



Po 



■go. 



d-v = 



(2.6) 
(2.7) 



In the first equation p a = (pi + p 2 )/2 and v is the kinematic viscosity. The quantity 
—<f>d(8T '1 8(f) I ' Po is the coupling term that accounts for capillary forces. It is easy to 
verify that it can be rewritten as —A (d 2 cf)d(f)) / p Q plus a gradient term which can be 
absorbed into the pressure term. Finally, p'g a / p a is the buoyancy contribution, p' being 
the deviation of the actual density, p, from the mean density p : 

p = p- Po ■ 

The buoyancy contribution can be rewritten in terms of p\ , P2 and <fi as: 

P - Po 



P 

—9o 

Po 



Po 



-g a 



pi 



P2 



Po 



Po 

= - A(j>g a 

where A = (p2 — Pi)/(p2 + Pi) is the Atwood number. 



■g a 



(2.8) 



2.3. Energetics 

Let us define the kinetic energy (per unit volume), Ek, and the potential energy (per 
unit volume), Ep, for our system ruled by Eqs. (|2.5| . (|2.6| and (|2.7| . 
By definition of potential energy, we have: 



Er 



dx dy p 2 gy- 
Pi)g 



1 



~PoAg(y 4>) 



dxdyp x gy- 



1 



E° P 



(2.9) 



being the total volume occupied by the fluids and brackets, (•••), denote spatial aver- 
ages. In Eq. (|2.9[) the constant E° p is chosen such to set the potential energy to zero for 
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vanishing Atwood number. 

In a similar way, one can define the kinetic energy per unit volume as: 



= 92 



ax dy p 2 — — + K / / dxdypx — — 



2 2 n 

,2 n 



~2> + 



T 



2 9 

P°(y) - Po-4(0 y) 



(2.10) 



From Eqs. (|2.5p and (|2.6[) we immediately realize that such equations are left invariant 
under the simultaneous transformation g — > —g, <f> — > —(f). As a consequence, ((j)v 2 /2) =0 
and the resulting kinetic energy simply reads: 

„,2 



Ek=Po( Y ) ■ 
By defining Ejr = T/Q the total energy of the two-fluid system is 

E = Ep + Ek + Ejr 



(2.11) 



The equation for Ek is obtained by multiplying Eq. (|2.6|) by /9 ^a and then taking spatial 
average. We easily get: 

».2 



dBjf/dt - podti-j) = -p v{{d a vY) + P oAg{v4>) - A(v a (5 Q 0) (d 2 cf)) 



Let us now take Eq. (]2 . 5[) . multiply it by y, and take the average: 



d t {yct>) + (yd y (#)) = jA(yd 2 -d 2 cf> + 



= 



(2.12) 



(2.13) 



by translational invariance and Leibniz rule. We thus have: 

dEp/dt = -d t (p Ag(y<fi)) = - p Ag{v(, 

where we have used the fact that (yd y (vcj))) = —((d y y)v(f)) - 
variation is 

ST < 



(2.14) 

(vcf)). The free-energy 



dtT= 1 1 ^Tt dxdy = 



ST 



i.e., 



-7( 
-7( 
-7< 

d t E. 



■'•do :■ ( ^ 

2 



dx dy 



ST 



^ ~ I / 77 Vidi<j) dx dy — 



ST 



>fi- A 



(-9 2 0) + 



dx dy 



)Sl-A{e afi (d a <l>)(ap<l>))Sl 



-7( 



9 a 



ST 



- A(e Q/3 {dp4>)) 



(2.15) 



(2.16) 



where we have introduced the strain tensor e a p = (d a vp + dpv a ) /2 and assumed bound- 
ary conditions suitable to justify integrations by parts. 
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The energy balance takes then the form: 



dt{E K + E P + Ejr) = - Po v{{d a vf) - 7 < 



1 2 



(2.17) 



The global system in thus intimately dissipative, even for a vanishing kinetic viscosity. 
It is worth emphasizing the cancellation of A(e Q| a(d a (/>)(9/30)) by the kinetic and the 
free-energy contributions, due to exchanges between the velocity field and the interface. 

2.4. Dispersion relation for the phase-field model 

The aim of this section is to show that the well-known dispersion relation for gravity- 
capillary waves (jChandrasekhar 1 119611 ) can be easily obtained within the phase-field 
formalism. To do that, let us concentrate our attention on a two-dimensional problem 
and indicate by y the gravity direction. Moreover, we will assume heavier fluid to be 
placed below the lighter one, in a way to have a stable situation. For a given perturbation 
imposed to the interface, the problem is to determine how the perturbation evolves in 
time. 

Denoting by h(x,t) a small perturbation imposed to a planar interface, we can rewrite 

<f) as: 

y - h(x,ty 



= f 



(2.18) 



where h can be larger than e, yet it has to be smaller than the scale of variation of h 
(small amplitudes). 

Locally, the interface is in equilibrium, i.e.: 



f" = V'{f) , 

where V{4>) = {4> 2 — l) 2 /4e 2 . In this limit we have: 



/t = -A 



a 2 / 

dx 2 



A 
e 



,d 2 h f" fdh x J 

dx 



(2.19) 



(2.20) 



dx 2 e 

Linearizing Eq. (|2.6p for small interface velocity we have, neglecting the viscous term: 

p dtv = -d y p - 4>dyH - Agporf 1 ■ ( 2 -21) 
The integration in the vertical direction interpreted in the principle value sense 



q y := lim 

LJoo 



Pod t q y := lim 

L\oo 



v dy 
L 

L 



yields: 



//"^ - -ft 
Ox z e 



d 2 h 



dh 

dx 



- Agp I fdy 



Pod t q y = °-q^2 - 2Agp h , 



having used the relations J{f') 2 dy = 2 v / 2/3, / ff'dy = and 

lim / fdy = +2h . 



(2.22) 
(2.23) 

(2.24) 

(2.25) 



The height variation of the interface has to match the vertical fluid velocity, thus giving: 



dth = v(x, h(x, t), t) 



(2.26) 
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The last step is to relate the velocity at the interface with the integral q y . This is done 
by restricting to potential flows: 

v = di\) d 2 i(i = . (2.27) 
For y > 0, denoting with """ the Fourier Transform, we have: 

/>oo 

i/j(x,y,t)= e~ ky+lkx ip(k,t)dk + c.c. (2.28) 



v(x, y,t) = - / ke- ky+lkx tp(k, t)dk + c.c. (2.29) 
Jo 

poo 

q y (x,t) = -2 e lkx 4>(k,t)dk + c.c. (2.30) 
Jo 

/>oo 

v {mt) = - ke lkx ^{k, t)dk + c.c. (2.31) 



Therefore: 



fi (i»t) = Mt , (2.32) 



so that in k— space we have: 



kq y 



d t h = J f Pod t q v = {-(rk z -2Ag Po )h . (2.33) 
From these two equations we immediately get: 

dfh + uj 2 h = , (2.34) 

with: 

uj 2 (k) = +Agk + — k 3 (2.35) 
2p Q 

that is the expected dispersion relation (jChandrasekhar 1 119611 ). For the stable config- 
uration we have, for all values of a: Agk + a / (2p ) k 3 > 0, i.e. any initially imposed 
perturbation will not grow indefinitely. 
From Eq. (|2.34|) and the initial condition: 

d t h(k,t) = att = , (2.36) 

we immediately have: 

h(k,t) = h(k,0) cos (cut) (2.37) 
and the velocity at the interface reads: 

v y nt (k,t) = -h(k,0)u sin (ut) . (2.38) 

Assuming an initial perturbation of the form h(x, 0) = ho cos (kx), from Eqs. (|2.3ip and 
(|2~38)) we obtain: 

$(k, t) = Q)u sin(cji) , (2.39) 

k 

and the velocity components, for y > 0, read: 

w T (x, y, t) = v(x, y,t) = - cos (kx)e~ ky h oj sin (cot) (2.40) 
u~(x, y, t) = u(x,y, t) = sin (kx)e~ ky h ui sin (uit) , (2-41) 
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where we used the relation ho = 2h(k, 0). 

For y < 0, in a similar way we obtain the velocity field components: 

v l (x,y,t) = v(x,y,t) = - cos (kx)e +ky h u sin (uxt) (2.42) 

u l (x,y,t) = u(x,y,t) — - sin (kx)e +ky h uj sin (ut) . (2.43) 

When in the initial configuration the heavier fluid placed above the lighter one, the 
dispersion relation (|2.35p trasforms in: 

u 2 (k) = -Agk + ^-P , (2.44) 

which is readly obtained by flipping the sign of g. For a < a c = 2p Q /(Agk 2 ) surface 
tension is not able to contrast gravity-induced vertical motion with the final result that 
amplitude perturbations grows exponentially: the flow is unstable. More precisely, from 
relation (|2.44j) and for a < a c we have: 



Lj(k) = \ -Agk + — k 3 = ia(k) , (2.45) 
V 2 p 

and Eqs. (f2T40|) - (|2~43|) transform in: 

v^(x,y,t) = v(x,y,t) = cos (kx)e~ kv hoasmh (at) (2-46) 
u^(x, y, t) = u(x, y, t) — — sin (kx)e~ ky hoa sinh (at) , (2.47) 



for y > 0, and: 



for y < 0. 



v l (x,y,t) ee v(x,y,t) = cos (kx)e +ky h a sinh (at) (2.48) 
u^(x,y,t) ee u(x,y,t) = sin (kx)e +ky h a sinh (at) , (2.49) 



3. Numerical investigation 

In this section we report results we have obtained exploiting direct numerical simula- 
tions (DNS) of the phase-field model for the Rayleigh-Taylor problem described in the 
preceeding sections. Our attention will be focused both on the linear phase of the per- 
turbation evolution and on the weakly nonlinear regime governed by plumes, for A 1. 
In the present study we will consider initial perturbations imposed to the interface vary- 
ing along one of the horizontal directions, say the cc-axis, and invariant along the other 
horizontal direction, say the z-axis. The perturbation is thus intimately two-dimensional 
a fact that allows us to solve the original Navier-Stokes equations coupled to the phase 
field in two dimensions. This clearly permits to obtain high accuracy and thus to properly 
test the phase-field approach against known results for both the linear and the nonlinear 
evolution stage. 

For a two-dimensional flow it is convenient to introduce the vorticity field u> [lo — (d x v) z ] 
and study the equations 



\L0 - 



v ■ duj = 



-vd 2 u) 



-—dx (d 2 (j) deb) - A (d(b) x g 

Po 



(3.1) 



d t (p + v ■ 8(f) = jd 2 n = 7 Ad 2 



(3.2) 
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In order to efficiently and accurately solve those equations we exploit a pseudospectral 
method |. Accordingly, periodic boundary conditions have to be assume d along the two 



direct i ons. For the horizon tal direction it is a natural choice (see e.g. ICabot fe Cook 



(|2006F l: lLiu fe Shen I l)2003l) ) while along the vertical one this choice deserves some com- 



ments. As initial condition we started from the hyperbolic-tangent profile, Eq. (|2.3p . for 
<p with the interface placed in the middle of the domain. The fact that we have periodic 
boundary conditions along y simply means that far from the middle of the domain the 
hyperbolic-tangent profile has to be distorted in order to satisfy periodic boundary con- 
ditions. However, both in the linear and in the weakly nonlinear regimes the amplitude 
of the interface perturbation is always much smaller than the vertical size of the box, 
so that the actual choice of boundary conditions at the top and bottom can be safely 
neglected. 

Such a strategy has been already exploited for the miscible case bv lCelani. Mazzino fe Vozella 
(|2006l) . 



The box has a horizontal to vertical aspect ratio L x /L y = 1 for the linear analysis stage 
and L x /L y = 1/2 for the weakly nonlinear evolution. In the latter case we take a smaller 
aspect ratio owing to the fact that the perturbation can reach a higher amplitude (with 
respect to case of the linear analysis) . 

In both cases the resolution is 1024 x 1024 collocation points. We need such a high reso- 
lution (despite the fact that we focus on a linear and weakly nonlinear study) in order to 
have a well described interface separating the two phases. In our simulations the mixing 
width (~ 4 e) is 6 mesh points. 

The time evolution is implemented by a standard second-order Runge-Kutta scheme. 
The physically relevant parameters in the present problem are the kinematic viscosity v, 
the buoyancy intensity Ag and the surface tension a. Both Ag and v will be varied in our 
study, while a will be kept fixed to a fixed value (see below) . The surface tension is related 
to the ratio A/e with e (and thus A) sufficiently small in order to have a finite value for 
the surface tension and, at the same time, to reproduce the correct sharp- interface limit. 
Finally, the parameter 7 appearing in the relaxation term in Eq. (|3.2p must satisfy the 
requirement that 7A be small, so as to enforce 'istantaneous' local equilibrium between 
flow and interface. Here we used the value (model units) 7A = 10~ 8 . 

All simulations presented here start from an initial condition corresponding to an 
equilibrium configuration: velocity identically zero and hyperbolic tangent profile for the 
phase field </>, expressed by the relation of the form: tanh ((y — h(x 1 t — 0))/c) with 

h{x. t = 0) = ho sin (fc x) 

For a given fc we choose the initial amplitude ho in a way that ho / X (where A = 27r/fc) 
is sufficiently small to fall in the linear phase (i.e. ho / A <C 1 ) and ho is sufficiently large 
for the wave disturbance to see an almost infinitesimal mixing width (i.e. ho / e 1 ). 
Specific numerical values are reported in the next sections. 

3.1. Linear instability for negligible viscosity 

The aim of this section is to verify the growth-rate (12.45P which holds in the linear phase 
when the viscosity is negligible. 

In order to do so, we take a small value of v {y = 10 -5 in the model units) and vary 
fc (up to fc c = {2Agp / a) 1 / 2 , the critical wave-number separating unstable from stable 
wave-modes) and Ag and take a fixed value of a. The ratio ho / A = 0.06 while ho / e 
ranges from ~ 10 to ~ 40 in the range of fc considered. 

The behavior of the square growth-rate a 2 is shown in dimensionless form in Fig. [2] as 
a function of fc for three different values of fc c (obtained by varying Ag) and in Fig. [3] 
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Figure 2. The square growth-rate a 2 (see Eq. (I2.45p ) for three different values of Ag corre- 
sponding to three different values of the critical wave number k c = (2Agp /o') 1 ^ 2 : k c — 3.4 
(solid circle), k c = 4.7 (solid triangle) and k c — 5.7 (solid rhombus). The dashed line is the 
linear-theory prediction expressed by the relation (|2.45fl . 
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Figure 3. The square growth-rate a 2 for k = 1 (solid circles), k = 2 (solid triangles) and k = 3 
(solid rhombus), all smaller than k c , for six different values of Ag ranging from 0.11 to 0.61. 
The dashed line corresponds to the linear-theory prediction. 



by varying Ag for three different values of k < k c . In both figures, symbols refer to the 
numerical results and the dashed line is the theoretical expectation given by (|2.45p . 
The numerical data in Figs. [5] and [3] have been obtained via best-fit of (v 2 ), the spatial 
average of v 2 as a function of time. The latter average is computed over a horizontal 
strip containing the interface (placed in the middle of the computational domain) and 
having an extension of a y above and below the interface. This has been done to avoid 
spurious contaminations coming from the upper and lower domain regions affected by 
the boundary conditions. In formulae: 

( W 2 ) = -L-L f dy [ dx(v l ) 2 + ^—^- [ " dy [ dx(v^) 2 

Za y Li x J- ay Jo ^ a y l^x Jo Jo 

= — — [-e- 2k a y + l] a 2 h 2 sinh 2 (a t) , (3.3) 

where we used the expression (12.46[) and (|2.48[) for and ir, respectively. 
The best fit has been done with a as unique free parameter and its high accuracy can be 
verified in Fig. 0] where we show the time evolution of (v 2 ) for k c — 4.7 (solid triangles 
in Fig. [5]) and for four values of k smaller than k c . At ta > 1.5 nonlinear effects start 
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Figure 4. Time behavior of (v 2 ) for k c = 4.7 (in Fig. [2] corresponding to the solid triangle) and 
for four values of k < k c . (a) k = 1, (b) fc = 2, (c) fc = 3 and (d) k — 4. The numerical results 
(symbols) are compared with the corresponding best fit expressions (see the text for details). 
In the insets the interface perturbation, h(x,t), is plotted at ta = 1.5 revealing a very accurate 
linear analysis prediction. 



to enter into play giving rise to corrections to the linear analysis (see Sec. 13. 4p . Up to 
that time, linear theory is very accurate as one can also realize by looking at the insets 
of Fig. [2] where the sinusoidal form of h(x,t) is reported for ta = 1.5. 

3.2. Linear instability for finite viscosity 

The aim of this section is to investigate numerically how the growth-rate, a, is modified 
by viscosity. As discussed in Appendix El both an upper and a lower bound for the per- 
turbation growth-rate are known (see Eqs. (|A lj) and (|A 2j) ) and we want to assess how 
the actual growth-rates compare with those. 

For such purpose, we choose a surface tension, a, and Ag in such a way to obtain in- 
stability for few (unstable) wavenumbers. Our choice was k c = 5.7 (see Sec. 13- 1[) thus 
corresponding to 5 unstable wavenumbers. 

As far as the initial perturbation is concerned, we report here the case corresponding 
to k = 1. Initial perturbations with a larger wavenumber simply need an initial smaller 
amplitude (and eventually a larger numerical resolution) in order to satisfy fto > e and 
ho -C A. Here, we have ho/X = 0.03 and hn/e ~ 20. Such ratios turned out to be suf- 
ficiently 'asymptotic' to produce accurate results. The effect of viscosity is studied by 
considering twelve values of viscosity in the range 10~ 5 < v < 5 10~ 2 (model units). 
The results of our simulations are summarized in Fig. [5] where the behavior of the square 
perturbation growth-rate, a 2 , is shown as a function of viscosity. The numerical predic- 
tions have been compared with the available theoretical bounds (dashed lines). 
Note that the numerical points are always in between the two bounds and also how the 
relative differences between the upper bound and the nu merical values are < 11 %. This 
latter fact is compatible, for example, with the results of Menikoff et al. I (|l977l h 
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Figure 5. Behavior of the dimensionless perturbation growth-rate, a v , for k = 1 and Ag cor- 
responding to k c = 5.7. Dotted lines correspond to upper and lower bounds for the growth-rate 
(see Eqs. (|X2|) and JXlJ). The arrow selects a value of the viscosity for which the time evolution 
of (v 2 ) is reported in the inset. The continuous line is the best fit slope (see text). 

The value of the growth-rates have been obtained via best of (v 2 ) (see Eq. (|3.3p ). Un- 
like what we did in previous section, here we perform the fit within the exponential 
region. The reason is that the non-asymptotic form of the perturbation time-evolution is 
unknown in the present case. 

The fit accuracy can be appreciated in the inset of Fig. [5] where the temporal evo- 
lution of the pertubation for v = 0.3 (model units) is shown together with the best fit 
slope (dashed line) from which a v is determined. Error bars, estimated by looking at the 
fit sensitivity by varying the length of the fit interval, are of the order of the symbol sizes. 

3.3. Stable configuration: gravity- capillary waves 

The performance of the phase-field approach in the unstable regime predicted by linear 
theory both in the presence and in the absence of viscosity proved to be very good. 
As discussed in Sec. 12 .41 for sufficiently large surface tensions and/or sufficiently small 
differences between fluids density, a perturbation initially imposed to the fluid interface 
may maintain its initial amplitude giving rise to the dispersion relation (|2.35j) . The waves 
resulting from the balance between gravity and surface tension are known as gravity- 
capillary waves. Our aim here is to verify their dispersion relation. 

To do that, we have fixed the parameters to obtain a critical wavenumber of order one. 
For Ag — 0.008 (model units) and the same a as in the unstable case, one has k c = 0.9. 
The first accessible wavenumber is thus stable and should evolve in time according to 
(|2.35p . However, the geometrical/computational configuration used in the unstable case 
did not produce sufficiently accurate results. In particular, using the same domain aspect 
ratio L x /L y — 1 and the same ratio between perturbation amplitude and perturbation 
wave-length we found a dynamics too dissipative with respect to what is expected. In 
the absence of viscosity, dissipation arises in the phase field formulation due to the 
sole contribution proportional to 7 in Eq. (|2.17j) . The latter parameter has been taken 
sufficiently small to ensure a negligible effects inside a period of oscillation. The specific 
value was 7 = 6.25 x 10~ 5 . To avoid spurious dissipation, as that induced by nonlinear 
effects, we reduced the amplitude of the initial perturbation with respect to the unstable 
case. Also, we increased the size of the periodicity box along the gravitational direction 
in a way to reduce possible spurious contribution arising from the upper/lower part of 
the computational domain where instabilities, not present in the unstable case, might 
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Figure 6. Time behavior of the perturbation maximum, n(t), for k = 1 and ho/X = 0.012. The 
critical wave number is k c = 0.9. Numerical results (symbols) are compared with the prediction 
from linear theory (see Eq. (|2,37[) ). 



now develop. The above choice on the amplitude of the inital perturbation implies a 
consequent reduction of e. The following set of parameters have been used: e = 0.008 , 
L x /L y = 1/4 and a resolution Nx x Ny of 256 x 4096. For an initial perturbation on 
k = 1, its initial amplitude h has been chosen to have ho/X — 0.012 and h /e ~ 10. The 
behavior of the maximum, n{t), of the initial perturbation is shown as a function of time 
in Fig. [6] The continuous line is relative to a sinusoidal with pulsation uj obtained from 
(|2 . 35|) . The agreement between theory and numerics is satisfactory both for the amplitude 
and for the pulsation. Note the small reduction of rj(t), in one oscillation period: only 1 
grid box over 4096. 

3.4. Weakly non-linear stage 

In this section we investigate the early stages of the nonlinear dynamics. We focus on 
the rising/falling velocity of plumes in the limit of small Atwood numbers when spikes 
and bubbles are known to coincide. The theoretical prediction for the terminal velocity 
is reported in Appendix [B] Our aim here is both to verify the existence of a regime 
characterized by a costant 'terminal' velocity and, secondly, to compare the prediction 
(|B lj) for such terminal velocity with our numerical data. 

The physical parameters are chosen to magnify the effect of the surface tension on the 
terminal velocity. This happens when the wavenumber k of the initial pertubation (still 
supposed unimodal) is slightly below k c . Here we choose Ag and a such that k c = 
4.004 and thus look at the dynamics associated to the wavenumber k = 4. The initial 
perturbation has an amplitude ho/X — 0.06; the initial dynamics is thus linear. Although 
we are interested to investigate the case of zero viscosity, in order to prevent numerical 
instabilities we add a small viscosity v — 2x 10~ 5 (model units). In Fig.[7]the perturbation 
amplitude is shown as a function of time: symbols correspond to our numerical data and 
the dashed line is the prediction (IB 1|) . A good agreement is found between numerics 
and theory in the range 1.2 < tU/X < 1.8. At larger times, neighboring plumes start to 
interact and the arguments leading to (|B lj) do not apply any longer. In Fig. [8] we show 
some snapshots of the evolution of the two fluids. Figures are equally spaced in time in 
the interval 1.2 < tU/X < 1.8. Black corresponds to (f> = —1; white to </) = 1. Their 
shape is similar to that experimentally observed. Note the aforementioned spike/bubble 
sy mmetry corresponding to the up-down symmetry of our original evolution equations. 



sy mmetry corr esponding to tne up -down ; 
bv lWaddell. Nieserhaus fc Jacobs I (|200lh . 
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Figure 8. Two-color snapshots of the phase field. Black (white) corresponds to (f> — — 1 
!> = 1). Frames are equally spaced in time in the interval 1.2 < tU/X < 1.8 (see also Fig. [7}. 



4. Conclusions and perspectives 

In this paper we showed that the phase-field model provides a valuable numerical in- 
strument for the study of immiscible, convective hydrodynamics. As a testground for this 
model, we have considered the Rayleigh-Taylor instability. Numerical results compare 
very well with known analytical results both for the linearly stable and unstable case, 
and for the weakly nonlinear stages of the latter. 

All these results are very encouraging in view of the next important step that is the the 
numerical simulation of immiscible RT turbulence. There, the interplay of all the funda- 
mental mechanisms that we have illustrated here (instabilities and wave propagation) is 
expected to give rise to a small-scale emulsion-like phase dominated by gravity-capillary 
waves and by a large-scale hydrodynamic range of scales where classical Kolmogorov tur- 
bulence should appear. This theoretical suggestion still awaits numerical confirmation, 
and the phase-field model provides the appropriate method to pursue this goal. 

We acknowledge useful discussions with Hekki Haario. AM and LV have been partially 
supported by PRIN 2005 project n. 2005027808 and by CINFAI consortium (AM). LV 
acknowledges support from From Discrete to Continuous models for Multiphase Flows 
TEKES project n. 40289/05. 
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Appendix A. Bounds for the perturbation growth-rate in the 
presence of viscosity 

The effect of viscosity is to reduce the perturbation growth-rate. However it does not 
remove the instabilities. Analytically, it is more difficult to consi der the effect of viscosit y 
with respect to surface tension (see Eq. (115) at page 443 of IChandrasekhar I Il96ll ). 
Nonetheless, it is possible to determine a lower and an upper bound to the growth-r ate 
a v . These bounds are the solutions to the following equations (jMenikoff et al. II1977T ): 



at + 2vk 2 ai + (v z h A - — )kat - fi/V + — )vk A a u - (^ 4 fc b - — )k A = (A 1) 

K k k £ 

al + 2vk 2 a u - a 2 = .(A 2) 

where a is the growth-rate in the inviscid case (see Eq. Q2.45p ), The solution of Eq. (|A 2[) 
is: 

a u = —k 2 v 



Vk^7 2 



(A3) 



while only a numerical solution is available for Eq. JXTJ. 

The goodness of those upper and lower bounds are numerically investigated in Sec. 
by means of the phase-field model. 



Appendix B. Models for the terminal bubbles/spike velocities in the 
weekly nonlinear regime 

Substantial deviations from the linear theory ar e obseryed w hen the perturbation am- 
plitude reaches a size of the order of 0.1 A - 0.4 A (jSharo 111984^ 



In that case the perturbation evolution is nonlinear. Then the disturbance grows non- 
linearly and the interface starts to deform. Indeed, at least for finite values of A, the 
interface can be divided into spikes corresponding to the regions where the heavier fluid 
penetrates into the lighter one, and bubbles associated to those regions where lighter 
fluid rises in the heavier one. The roll- up of vortices produces a mushroom- type shape 
for bubbles and spikes (see, for instance, Waddell. Nieserhaus fc Jacobs 112001 ). When the 
fluid densities are similar (corresponding to our case A <C 1) spikes and bubbles coincide 
and approach a constant and equal velocity. In both cases, the exponential growth of 
the velocity perturbation amplitude c haracterizing the linear phase of the evolution is 
replaced by a linear-in-time behavior (jWaddell. Nieserhaus fc Jacobs 1120011). Two mod- 



els are available to d escribe this stage: the drag-buoyancy model (lAlqn et al. 1995) and 
the "Layzer model" (jLavzer Ill955l lGoncharov1l2003t lYoung fc Ham 1120061 ). The former 
model describes bubble and spike motion by balancing the buoyancy and drag forces and 
it assumes that this velocities reach a constant values for sufficiently long times. The lat- 
ter model uses an expansion of the perturbation amplitudes and conservation equations 
near the tip of bubbles an d spikes. Thi s approach has been first applied to the fluid- 
vacu um interface (A = 1) (Layzer. 19551 ) and then extended to arbitrar y Atwood num- 



bers ( Goncharov 20031 ) and to include the surface tension contribution ( Young fc Ham 
2006). According to the latter study, in our case (bidimcnsional flow, immiscible fluids 



and smal l Atwood number) on e expects that the terminal bubble and spike velocity be 



equal to (|Young fc Haml l2006): 



U(t 



oo 



-A- - - 

3 k 9 p 2 



Pi 



(Bl) 



This expectation is numerically tested, in Sec. 13.41 by exploiting the phase-field method. 
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